;******************************************
; plot regions with climate regime change 
;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;********************************************
 begin
 
 dirPath = ".../NFood_replication/data/" ; change to local path of your directory

 nmon = 12;
 syr = 2006;
 fyr = 2100;
 csyr = "" + syr
 cfyr = "" + fyr
 nyr = fyr - syr + 1
 dmis = 9.96921e+36

 case = (/"RCP85","RCP45","RCP85_SAI","RCP85_MSB","RCP85_CCT"/);
 ncase = dimsizes(case);

;read dimension
 fdim = addfile(dirPath+"/land_cover/RCP85.1deg.landarea.nc", "r");
 lat = fdim->lat
 lon = fdim->lon
 nlat = dimsizes(lat);
 nlon = dimsizes(lon);
 dims = (/nlat,nlon/);
;print(lon)

;read data 
 dd2010 = new((/ncase,nlat,nlon/),"float"); Control period 2006-2015 for RCP85, RCP45
 ;dd2010@_FillValue = default_fillvalue("float")
 dd2050 = new((/ncase,nlat,nlon/),"float");
 dd2100 = new((/ncase,nlat,nlon/),"float");
 pc2050 = new((/ncase,nlat,nlon/),"float");
 pc2100 = new((/ncase,nlat,nlon/),"float");
 do kk=0,ncase-1
   if (kk .lt. 2) then; ;RCP45 and RCP85
     fi := addfile(dirPath+"/climate_data/"+case(kk)+".clm2.h0.2006-2100.climate.nc", "r");
     dd := fi->RAIN + fi->SNOW ;precipitation
     dd0 := dd(0:119,:,:)   ;CNTL time period: 2006-2015
     dd1 := dd(480:587,:,:) ;future time slice 1: 2046-2055
     ;dd2 := dd(1020:1127,:,:) ;future time slice 2: 2091-2099
     dd2 := dd(828:1127,:,:) ;future time slice 2: 2075-2099
   else
     fi := addfile(dirPath+"/climate_data/"+case(kk)+".clm2.h0.2020-2100.climate.nc", "r"); 
     dd := fi->RAIN + fi->SNOW
     dd1 := dd(312:419,:,:) ;future time slice 1: 2046-2055 
     ;dd2 := dd(852:959,:,:) ;future time slice 2: 2091-2099
     dd2 := dd(660:959,:,:) ;future time slice 2: 2075-2099

   end if
   dd2050(kk,:,:) = tofloat(dim_avg_n_Wrap(dd1,0)) ;10-yr avg along the first dim
   dd2100(kk,:,:) = tofloat(dim_avg_n_Wrap(dd2,0)) ; 
   pc2050(kk,:,:) = dd2050(kk,:,:)-dd2050(0,:,:) ; change relative to RCP45 
   pc2100(kk,:,:) = dd2100(kk,:,:)-dd2100(0,:,:) 
 end do
 

;******************

 pptch = pc2100* 86400*365 ;change in ppt, mm/s to mm/year
 ;printVarSummary(pptch) 
 printMinMax(pptch,0) ; range 
 ;min=-1567  max=2300
 ;print("mean="+avg(pptch))
 opt = True
 opt@PrintStat = True
 cstat = stat_dispersion(pptch, opt) ;detailed statistics
 ;set color range between lower/upper 5% and min/max
 ;set color range to -1000~1000 (mm) 

 pptch!0="case"
 pptch!1="lat"
 pptch&lat = lat
 pptch!2="lon"
 pptch&lon = lon
 pptch1 = lonFlip(pptch); 0~360 to -180~180
; copy_VarCoords(dd, ww);


;plot
 imagePath = dirPath+"/../figures/supplementary"
 imagename = imagePath+"/Fig.S2c"
 wks = gsn_open_wks("pdf",imagename)
 ;gsn_define_colormap(wks, "ncl_default")             ;choose a color palette

 plot1 = new(ncase,graphic)

 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False
 res@lbLabelBarOn         = False                 ;turn off individual label bars 

 res@cnLineLabelsOn       = False                 ; turn off line labels
 res@cnFillOn             = True                  ; turn on color fill
 ;res@cnFillPalette        = cmap                 ; cnFillPalette spans the color map by default
 ;res@gsnSpreadColors      = True                 ; this based on an old color scheme that was deprecated
 res@lbLabelAutoStride    = True                  ; optimal labels
 res@cnLinesOn            = False

 res@gsnStringFontHeightF  = 0.020
 res@tmXBLabelFontHeightF  = 0.020
 res@tmYLLabelFontHeightF  = 0.020

; res@cnFillMode = "RasterFill"

 res@mpCenterLonF = 0.
 res@mpMinLatF = -60.
 res@mpMaxLatF =  60.

 res@mpShapeMode  = "FreeAspect"
 res@vpWidthF      = 0.8   
 res@vpHeightF     = 0.4   

 res@tmXBMode = "Explicit"
 res@tmXBValues = (/-180., -120., -60., 0., 60., 120., 180./)
 res@tmXBLabels = (/"180W","120W","60W", "0", "60E","120E","180E"/)
 res@tmYLMode = "Explicit"
 res@tmYLValues = (/-60., -30., 0., 30., 60., 90./)
 res@tmYLLabels = (/"60S","30S","0","30N","60N","90N"/)


; to have a common label bar, all plots should be set to the same interval
; b/c the label bar is drawn by default from the interval of the first plot.
  res@cnLevelSelectionMode =  "ManualLevels"   
  res@cnMinLevelValF       = -200  ; mm
  res@cnMaxLevelValF       = 200   ;  
  res@cnLevelSpacingF      = 20 ;       

; change a default color map
  ;---Reading the default colormap into a two dimensional 254 x 4 (RGBA) array
  cmap = read_colormap_file("NEO_div_vegetation_a") 
  ;cmap = read_colormap_file("ncl_default")
  ;cmap = cmap(::-1,:) ; reverse color map blue for wet, red for dry
  ;select the middle cmap elements that correspond to 0 in the middle of value range (-0.5~0.5)
  cmap(120:135,:) = 1 ; use white for zeros; 0 for black, 1 for white
  ;cmap(0,:) = 1 ; use white for zeros; 0 for black, 1 for white
  res@cnMissingValFillColor = "white"  ; missing values (default gray)
  res@cnFillPalette         = cmap     ; use the modified color map; cnFillPalette spans the color map by default

; res@cnLevelSelectionMode = "ExplicitLevels"
; res@cnLevels    = (/-1., -0.75, -0.5, -0.25, -0.001, 0.001, 0.25, 0.5, 0.75, 1./) ; level to divide data into ranges
; ;res@cnFillColors = (/ 0, 54, 2, 3, 4, 5, 28, 20, 14, 0/)
  ;res@cnSpanFillPalette     = True    ; spans cnFillPalette to get color index values for cnFillColors in Explicit mode

 scenario = (/ "RCP45 - RCP85", "SAI - RCP85", "MSB - RCP85", "CCT - RCP85" /)
 do kk=1,ncase-1
   res@gsnLeftString      = scenario(kk-1)
   plot1(kk)              = gsn_csm_contour_map_ce(wks,pptch1(kk,:,:),res)
 end do


;panel
 panres1                  = True
 panres1@gsnFrame         = False
 panres1@gsnPanelBottom   = 0.05          ;add some pace in the bottom
 panres1@gsnPanelYWhiteSpacePercent = 5   ;add 5% white space between plots
 panres1@gsnPanelMainString = "Precipitation change relative to RCP85 (2075-2099)"     ; set panel title
 panres1@gsnPanelLabelBar = True ;add a common label to whole panel
 panres1@lbLabelFontHeightF = 0.011 ;set label tick font size
 gsn_panel(wks,plot1(1:4),(/2,2/),panres1)        ; now draw as on

; Draw a text string at the bottom below label bar
  txres               = True
  txres@txFontHeightF = 0.015
  gsn_text_ndc(wks,"~F8~D~F~Precipitation (mm year~S~-1~N~)",0.5,0.15,txres) 
  ;function codes: ~F8~D is for big delta; ~F~ or ~N~ means switch back to default font
  ;gsn_text_ndc(wks,"~F8~D~",0.5,0.1,txres) ;~F8~D~ is the special code for delta

 frame(wks)


end
